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The quantum nature of nuclei plays an important role in the accurate modelling of light atoms 
such as hydrogen, but it is often neglected in simulations due to the high computational overhead 
involved. It has recently been shown that zero-point energy effects can be included comparatively 
cheaply in simulations of harmonic and quasi-harmonic systems by augmenting classical molecular 
dynamics with a generalized Langevin equation (GLE). Here we describe how a similar approach 
can be used to accelerate the convergence of path integral (PI) molecular dynamics to the exact 
quantum mechanical result in more strongly anharmonic systems exhibiting both zero point energy 
and tunnelling effects. The resulting PI-GLE method is illustrated with applications to a double-well 
tunnelling problem and to liquid water. 



I. INTRODUCTION 

Atomistic computer simulations have become an im- 
portant complement to experimental measurements in 
studies of a wide variety of physical, chemical and biolog- 
ical systems. However, in order to reduce computational 
effort that is required for larger systems, a number of 
approximations are typically made in these simulations. 
One that is commonly employed is to assume that the 
atomic nuclei behave as classical particles, even when the 
interactions between them are obtained from an ah initio 
calculation. This is a reasonable assumption for heavy 
atoms at high temperatures. But for lighter atoms - and 
hydrogen in particular - significant deviations are to be 
expected from classical behavior even at room tempera- 
ture. 

When empirical interaction potentials are used to com- 
pute the forces acting on the nuclei, nuclear quantum 
effects are often implicitly accounted for, because force 
fields are typically parameterised so as to agree with 
experimental data when used in classical molecular dy- 
namics simulations. When ah initio methods are used 
to describe the interactions between the nuclei, there is 
no such parameterisation for nuclear quantum effects, 
and the error that results from assuming purely classical 
behavior of the nuclear motion is often comparable to 
that stemming from an approximate treatment of elec- 
tron correlation.ti] 

The conventional approach to including nuclear quan- 
tum effects exploits the isomorphism between the quan- 
tum mechanical partition function of the physical system 
and the classical partition function of an extended prob- 
lem consisting of a necklace (or ring polymer) of repli- 
cas of the system in which corresponding atoms are con- 
nected by harmonic springs.^ As the number of replicas 
(or beads of the necklace) is increased, this imaginary 
time path integraP (PI) method samples an ensemble 
that converges systematically to that of the quantum 



problem with distinguishable particles, at the expense 
of a computational effort that increases in proportion to 
the number of beads. 

Methods have been devised to reduce this effort by 
splitting the calculation of forces into an inexpensive, 
short-ranged contribution that is computed on every 
bead, and a long-range contribution that is evaluated on 
a contracted ring polymer with fewer beads .'^''S! These 
methods are straightforward to implement in simula- 
tions with empirical force fields. However they cannot 
presently be used in ah initio simulations, where the 
splitting of the forces into an inexpensive short-range 
contribution plus a more expensive long-range contribu- 
tion is significantly more difficult to arrange. Approx- 
imate quantum methods such as the Feynman-Hibbs'^-'^ 
approach are also difficult to use in ah initio simulations, 
because they require the computation of the Hessian. 

A more general approach to reducing the effort of path 
integral simulations is therefore called for, and there are 
certain indications in the literature that this might now 
be possible. In particular, it has recently been shown that 
colored noise, generalized Langevin equation (GLE) ther- 
mostats can be used not only to enhance the s amp ling 
of classical and path integral molecular dynamics,'^ but 
also to modify conventional molecular dynamics so as 
to include nuclear quantum effects in mildly anharmonic 
systems in which zero-point energy plays a significant 
role.^*^ This approach involves negligible overhead with 
respect to purely classical dynamics, provides very nat- 
urally the proton momentum distribution (which is rel- 
evant for comparison with inelastic neutron scattering 
experiments), and can be applied with the same ease to 
empirical and ah initio force fields. However, it does not 
provide a satisfactory description of more subtle quan- 
tum effects such as tunnelling, and its accuracy cannot 
be systematically improved. 

In this paper, we will show that it is possible to com- 
bine a GLE thermostat with path integral molecular dy- 
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namics (PIMD) in such a way as to exploit the best fea- 
tures of both techniques. In particular, by tuning the 
properties of the correlated noise in an appropriate GLE, 
we will show that the systematic convergence of PIMD to 
the exact quantum mechanical result can be greatly ac- 
celerated, leading to significant computational savings for 
a given level of accuracy. The resulting PI-I-GLE scheme 
provides an equally reliable description of zero point en- 
ergy and tunnelling effects, it is equally applicable to sim- 
ulations with empirical and ab initio force fields, and un- 
like the original single-bead "quantum thermostat" dis- 
cussed above it can be systematically improved simply 
by increasing the number of path integral beads. 

The outline of the paper is as follows. In Section II 
we briefly recall a few concepts from path integral and 
generalized Langevin equation methods, and describe our 
strategy to have them work in synergy. In Section III we 
present a systematic study of a one-dimensional quartic 
double well potential, discussing the effect of zero point 
energy and tunnelling on the convergence of our method. 
In Section IV we examine how the method performs for 
liquid water, and in Section V we draw our conclusions. 



II. COMBINING PATH INTEGRALS WITH 
THE GENERALIZED LANGEVIN EQUATION 

A. Path integral methods 

Let us first recall the basic principles of imaginary time 
path integral methods, so as to introduce our notation. 
We will only discuss the o ne-dim ensional case and refer 
the reader to the literatur^^EIIlHlI! for a more detailed 
discussion. 

Consider the Hamiltonian for a particle in a one- 
dimensional potential. 



H=-f 



(1) 



where p and q are the mass-scaled momentum and posi- 
tion operators and the potential V{q) is assumed to be 
such that the quantum mechanical partition function 



(2) 



is well defined. The path integral formalism avoids the 
solution of the Schrodinger equation for the Hamiltonian 
in Eq. ([!]) and allows one instead to sample configura- 
tions consistent with the quantum mechanical equilib- 
rium distribution by exploiting an isomorphism with an 
extended classical problem.^ Indeed a standard Trotter- 
producl)^ approximation to the Boltzmann operator in 
Eq. ([2]) yields the following extended phase space expres- 
sion for the partition function 



d^p / d^qe 



-Hp(p,q)/PfeBT 



, (3) 



where iJp (p,q) is the classical Hamiltonian of a ficti- 
tious ring polymer composed of P replicas of the system 
connected by harmonic springs 



ffp(p,q) 



p-i 

E 



1 



qj+i) 



(4) 

with ujp = PkBT/h and qp = qq. The error in this 
approximation is 0{P~^) and so vanishes to leave the 
exact quantum mechanical result in the limit as P — > 

The momenta are often integrated out of Eq. ([s]) 
to leave a purely configurational integral which forms 
the basis of the path integral Monte Carlo (PIMC) 
technique.!^ One can however retain the momenta, and 
describe the dynamical evolution of the ring polymer 
by mean s of H amilton's equations. This is the PIMD 
approachj ^^ l ^^ l which provides a particularly efficient way 
to sample the configuration space in situations such as 
molecular liquid simulations in which an effective PIMC 
calculation would require complicated collective moves. 
A fully converged PIMD calculation produces the exact 
thermodynamic and structural properties of the quan- 
tum mechanical system. Although it is something of 
an aside to the present work, in which wc shall be con- 
cerned exclusively with static equilibrium properties, it 
is also n ow we ll established that the centroid mole cular 
dynamic^i^E^ and ring polymer molecular dynamicJ^^^^ 
generalizations of PIMD can be used to provide quite 
reasonable estimates of dynamical properties such as dif- 
fusion coefficients and chemical reaction rates. 

The extension to three dimensions and to an arbi- 
trary number of interacting distinguishable particles is 
straightforward. But rather than describe this extension 
here, let us consider instead the simple case of a one- 
dimensional harmonic potential, V{q) — uPq'^ /2. The re- 
sult will be used in what follows. For this simple model, 
the integral in Eq. ([3| can be evaluated by transforming 
to the P normal modes {qk}^~^ of the ring polymer with 
frequencies 



Wfc = + 4a;p sin^ {kn / P) . 



(5) 



It is then easy to show that the equilibrium position dis- 
tribution of each bead of the ring polymer will be a Gaus- 
sian centred on q — Q, with a variance 



p-i 



p-i 



(6) 



fe=0 



fc=0 



which converges to the correct quantum mechanical ther- 
mal expectation value 



(q ) = — coth — 

/ 2u} 2kBT 



(7) 



in the limit as P — 00. 
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B. Generalized Langevin equations 

It has recently been shown that an appropriate gener- 
alized Langevin equation (GLE) thermostat can be used 
to sample configurations for a harmonic oscillator that 
are consistent with the quantum mechanical variance in 
Eq. ^ using just a single replica of the system, thereby 
avoiding the overhead of a P-bead PI simulation!^ 

The basic idea behind GLE thermostats is to con- 
struct a linear, Markovian stochastic differential equa- 
tion (SDE) in an extended momentum space which is 
coupled to the Hamiltonian dynamics in such a way that 
the equations of motion read 



q = p 



-V'{q) 




O-pp ^p \ I P 

a„ A 



b 

b„ B 



(8) 



where ^ is a vector of rt+l uncorrelated Gaussian random 
numbers with (^j {t)£,j (0)) = 6ij6{t). It can be shown 
that the resulting trajectories are statistically equivalent 
to those obtained from a non-Markovian Langevin equa- 
tion involving only p and 

q^p 

P = -V'{q)- I K{t~s)p(s)ds + at), 



where the friction kernel K{t) and the noise correlation 
function H{t) — (C(OC(0)) ^^re related by analytical ex- 
pressions to the matrices that appear in Eq. (|8|. These 
relationships, together with a more detailed discussion, 
can be found in Ref.*^ 

The only non- linear term in Eq. ([s]) is the force V'{q). 
If the potential is harmonic, the force depends linearly on 
q, and the whole set of equations describe an Ornstein- 
Uhlenbeck process,I23 which can be solved analytically 
to yield closed-form expressions for static and dynamic 
properties of the trajectory. Based on these expressions, 
one can iteratively refine the parameters Op, aj, a^, A, 
bp, bp , bp and B in Eq. (j8| subject to certain positivity 
constraints until the desired response of the thermostat 
is obtained.l^ 

When considering the generalisation of this strategy to 
a multi-dimensional harmonic problem, one realises that, 
because of the linear nature of the SDE and of the con- 
sequent rotational invariance, the dynamics will conform 
to the analytical predictions obtained from the Ornstein- 
Uhlenbeck process even if the thermostats are applied to 
Cartesian coordinates, without the need to diagonalise 
the Hessian and transform to normal modes. Moreover, 
it has been demonstrated in previous work that the ana- 
lytical properties of the thermostat obtained in the har- 
monic limit provide a meaningful prediction of its be- 
havior for anharmonic problems with a similar range of 
frequencies.'^ 



By employing this strategy it is possible to obtain a 
number of useful effects, such as efficient sampling of the 
canonical distribution in constant-temperature MD. In 
this application, a fluctuation-dissipation theorem must 
hold, which relates the friction kernel and the memory 
of the noise in Eq. (|9| through H{t) = kBTK{t). More 
generally, when this condition is not imposed, one ob- 
tains a non-equilibrium dynamics in which a frequency- 
dependent effective temperature is enforced. Either way, 
one can calculate the stationary covariance of the har- 
monic dynamics in the {q,p,s) extended phase space, 
and the mean squared fluctuations of the positions and 
momenta, which we will label Cqq{ui) and Cpp{uj), respec- 
tively. It is then possible to design a GLE dynamics 
which behaves as a "quantum thermostat" fi^l i_e_^ which 
enforces probability distributions of positions and mo- 
menta that are consistent with those expected for a quan- 
tum harmonic oscillator, 



,2 '^PP 



h Huj 

— coth , 

2a; 2kBT' 



(10) 



and does so over a wide range of frequencies. 

When applying this idea to a multi-dimensional sys- 
tem one faces the problem of zero-point energy leakage.^ 
Anharmonic coupling causes a flow of heat from high- 
frequency to low-frequency vibrations, and a departure 
from the desired behavior in Eq. (10). This problem 



can be mitigated to some extent by exploiting the flex- 
ibility of the GLE in Eq. ([8| to enhance the coupling 
strength of the thermostat and ensure that all vibra- 
tions are maintained at the correct effective temperature. 
This approach has been shown to give satisfactory results 
in a number of realistic, condensed-matter applications, 
ranging from the calculation of diamond-graphite coex- 
istence curve^^ll to the proton momentum distribution in 
hydrogen-storage materials.'23 

Being simple to implement and computationally inex- 
pensive, this "quantum thermostat" is an important step 
towards a routine treatment of nuclear quantum effects 
in molecular dynamics. It lacks however two desirable 
features; namely, the possibility of treating more subtle 
quantum effects such as tunnelling, and the possibility 
of increasing the accuracy in a systematic way. Since 
both these requirements are met by PI methods, one sus- 
pects that it might be possible to develop a synergistic 
approach which combines path integrals with the GLE 
thermostat so as to obtain accurate results without the 
effort of a fully-converged PI simulation. The implemen- 
tation of such a PI-f-GLE approach is the subject of the 
next section. 



C. A synergistic combination 

When a GLE thermostat is applied to a path integral 
molecular dynamics simulation, the internal frequencies 
of the ring polymer necklace will be present along with 
the physical vibrations of the system. It is therefore once 
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again instructive to consider a harmonic model with fre- 
quency w, for which both the ring polymer frequencies 
and the stochastic dynamics can be treated analytically. 

In order to construct a non-equilibrium Langevin 
dynamics that enforces the quantum mechanical equi- 
librium distribution corresponding to the frequency- 
dependent fluctuations in Eq. (10), one must allow in 



the path integral context for the fact that the average of 
is obtained from a sum over the ring polymer normal 
modes, 



p-i 



p-i 



(«%-:^E('^^)-^E^..M' (11) 

■ ' k=0 



k=0 



where the normal mode frequencies cuk are given in 
Eq. ([5]). Here the last equality holds if a GLE which 
results in the frequency-dependent position fluctuation 
Cqqiu)) has been applied separately to each bead of the 
ring polymer. 

It follows from this that one cannot simply tune the 
parameters in the GLE acting on each bead of the neck- 
lace so that Cqq{uj) is given by Eq. (10 1. Instead, the 



appropriate frequency-dependent distrTbution to be en- 
forced depends on the number of beads in the necklace, 
and can be obtained by solving 



1 

- Cqq{uJk) 
k=0 



— coth . 

2uj 2kBT 



(12) 
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FIG. 1. (a) Relative error in the estimates of gs{x) at dif- 
ferent iterations of the self-consistent procedure in Eq. ( |17[ l, 
with a — 1/8. The wiggles which appear after many itera- 
tions occur because each gp\x) is approximated as a spline 
interpolation before computing gp'^^^x). These wiggles can 
be systematically reduced by using a denser grid in x. (b) 
Converged gp(x) curves for dilferent bead numbers. Note 
that the curves with larger values of P are flat up to larger 
values of x. 



The first task is therefore to solve this functional equation 
for Cqq{ijj), recalling again that the frequencies ujk are 
related to the frequency w of the harmonic oscillator by 
Eq. ([5| ; the solution will be a universal function of to that 
is equally applicable to any harmonic oscillator, just as 
in the case of the original quantum thermostatP^ that 
enforces the condition on Cqq{uj) in Eq. (10 1. 

Before we describe our approach to solving Eq. (12 1, 
let us transform it into dimensionless form by defin- 
ing X = hiu/2kBT, h{x) — xcothx, and gp{x) — 
{kBT/P){2x/hfcqq{2xkBT/h). With these definitions, 
Eq. ([12]) becomes 



p 

E 

fe=0 



^ gpjxk) 

xllx^ 



h{x), 



where 



7-2 



-I- p2 sin^ 



~P' 



(13) 



(14) 



The solution to Eq. ( 13 ) is not unique, and one must pick 
a particular solution by means of appropriate boundary 
conditions. In particular, one would like to enforce a 
physical behavior on the solution, with no discontinuities 
and reasonable asymptotic forms in the a; — > and x — > 
00 limits. The tentative solution 



g'-^\x)^hix/P) 



(15) 



satisfies these requirements, and it is a good approxima- 
tion to gp{x) from several points of view. First of all, 
it is the exact solution in the one-bead case (which cor- 
responds to the bare quantum thermostat) and in the 
infinite bead limit, where it yields a constant effective 
temperature on all frequencies (which is correct, as in 
this case PIMD alone is sufficient to converge to the ap- 
propriate distribution). It also provides the appropriate 



large- a; limit for a smooth solution to Eq. (13), for arbi- 
trary P. 

In order to refine this tentative solution, one can cast 



Eq. ( 13 ) into a fixed-point iteration, by singling out the 



fc = term: 



p-i 



gp{x) = h{x) - 



k=l 



gpjxk) 
xllx^ ■ 



(16) 



We have found empirically that a self-consistent itera- 
tive procedure that yields an exact solution gp{x) to 
Eq. ( 13 ) to arbitrary precision can be obtained by stabil- 
ising Eq. (|16|) with a mixing strategy. 



gf\x)^h{xlP) 
g^p^'^^ {x) — a h{x) 



9p\xk) 

fe=i 



x^/x-^ 



{l^a)gf{x). 

(17) 
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FIG. 2. Comparison between the properties of the stochas- 
tic dynamics which result from the fitting strategy used in 
the present work (left column) and that used for the origi- 
nal quantum thermostalP^ (right column) , for simulations at 
T = 298 K. Note that by sacrificing the fit to the fluctuations 
of momenta [cpp{ij)] it is possible to obtain a perfect fit to 
Cqq{uj) and a better sampling efficiency kv('^) = l/[''"v('^)i^], 
where tv{oj) is the correlation time of the potential energy 
for a harmonic oscillator with frequency uu. In the case of 
PI-I-GLE it is less essential to enforce strong coupling to avoid 
zero point energy leakage, and it is therefore possible to avoid 
the overdamping which hinders the sampling in the case of the 
original quantum thermostat. This overdamping is clear from 
the bottom right-hand panel, which shows that K{lj) ^ cj at 
low frequencies where K{ijj) is the Fourier transform of the 
memory kernel in Eq. 



In particular, the mixing parameter a = 1/P was found 
to give a sufficiently fast and convergent iteration for all 
of the numbers of beads we tried. The convergence of 
gp\x) to gp{x) is shown in Fig. 1 for P = 8, along 
with examples of the resulting converged solutions for 
other values of P. Tabulated values of these solutions 
for P ranging from 1 to 128 can be downloaded from an 
on-line repositorj^^H. Once one has a converged gp{x), 
the frequency-dependent position fluctuation Cqq{u)) in 



Eq. ( 11 ) is given by 



c,,(w) - {PkBT/Lj^)gp{hLj/2kBT), 



(18) 



and the only remaining problem is to design a GLE that 
can be used to enforce this fluctuation on each ring poly- 
mer bead. 



D. Fitting and implementation 



The design of a GLE consistent with Cqq{uj) in Eq. ( 18 ) 
consists of optimising the matrices in Eq. ([8| until the 
desired response of the thermostat is obtained. As has 
been found for a variety of other applications of the 



GLE,lsHl0I26|the efficacy of the resulting PI-fGLE scheme 
will depend on the strategy by which the optimisation 
is performed, and on the additional criteria besides fit- 
ting Cqq{uj) that are used to define the merit function 
for the optimisation. For applications of the quantum 
thermostat to anharmonic problems the strength of the 
coupling between the thermostat and the Hamiltonian 
dynamics and the efficiency of the sampling can be just 
as important as the agreement between Cqq{ui) and the 
target function in Eq. (18). The possibility of improv- 
ing the accuracy systematically by increasing the number 
of beads makes zero-point energy leakage and other an- 
harmonic effects a lesser concern in the present context. 
However, if systematic convergence is to be achieved, it is 
advisable that the sampling efficiency does not differ dra- 
matically between the fits for different numbers of beads. 

In the case of the original quantum thermostat de- 
scribed in Ref.I^, we chose to constrain Cpp{uj) = 



Lo Cqq{uj) ss in Eq. (10 1. This had the advantage of pro- 
viding ready access to the quantum mechanical momen- 
tum distribution, which is a quantity of direct relevance 
to recent deep inelastic neutron scattering experiments. 
However, by constraining Cppicd) in this way, we found 
that we had to enforce a strongly overdamped regime 
at low frequencies in order to avoid zero point energy 
leakage in applications to multi-dimensional systems.^Sl 
In the present PI-I-GLE scheme, as in PIMD itself, the 
ring polymer momenta lose all physical significance as 
soon as there is more than one bead, and extracting the 
quantum mechanical momentum distribution requires a 
considerably more intricate calculation.!^ When it comes 
to fitting the GLE, it is therefore more natural to regard 
the momenta simply as a sampling device, and to focus 
exclusively on the fluctuations of configurations Cqq{co), 
as we have already done in our description of PI-I-GLE 
in Sec. III.C. 

As shown in Fig. 2, this leaves significantly more free- 
dom in the fit, even in the case of just one bead. The 
extra freedom allows us to reproduce the desired Cqq{u!) 
with a maximum discrepancy smaller than 0.5% and to 
achieve high sampling efficiency over a broad range of 
frequencies, yielding a sampling performan ce c ompara- 
ble to that of an "optimal sampling" GLE.I^ To give 
another example involving more beads, the self-diffusion 
coefficient for a model of liquid water (see Sec. IV) - 
which in this context can only be regarded as a mea- 
sure of the sampling efficiency for slow, collective mo- 
tion - is reduced by less than 50% in a well-converged (8 
bead) PI-I-GLE simulation compared to NVE dynamics, 
whereas the original (one bead) quantum thermostat de- 
creases the diffusion coefficient of the same model by a 
factor of 10. 

Having obtained a nearly-constant sampling efficiency 
is also beneficial to the transferability of the fitted param- 
eters, as discussed in Refs.!^. As a matter of fact, the 
very same parameterization has been adopted for both 
of the examples given below, which are as different as 
a one-dimensional quartic double well and liquid water. 
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With an appropriate scalingpl^ these GLE parameters 
can be safely adopted in all circumstances where the max- 
imum physical frequency present in the system is smaller 
than S^ksT/h. The GLE parameters we have used in 
the present study are available up to P = 16 and may be 
downloaded from an online repositorji^H. 



The details of how we actually optimised the GLE 
matrices for the present PI+GLE application are rather 
technical, and less important than the criteria employed 
for the optimisation that we have just described. A com- 
prehensive discussion of the optimisation of GLE ma- 
trices has recently been given elsewhere P and we used 
exactly the same techniques in the present study. The 
implementation of a GLE thermostat into a PIMD code 
has also been discussed in detail in a recent paperj^where 
it was shown to be essentially no more complicated than 
adding a thermostat to a classical molecular dynamics 
simulation. However, in the case of PI-I-GLE there are a 
few additional points that we do need to make. 

First of all, the dynamics in PI-I-GLE must be per- 
formed using physical bead masses, as opposed to the 
alternative masses that are sometimes used in path inte- 
gral schemes. Unless physical bead masses are used the 
frequencies of the necklace will be different from those 
in Eq. and Eq. (14 1 will have to be modified ac- 
cordingly; this will change the functional equation for 
gp{x) in Eq. (13) and the new equation will have to be 
solved from scratch. The equations of motion can be inte- 
grated efficiently with physical bead masses by perform- 
ing a normal mode transformation for the free ring poly- 
mer ev olutio n^ and/or employing a multiple time step 
scheme.l^IIll 

Another important observation is that, at variance 
with canonical sampling GLE schemes in which the free 
particle propagation of the GLE preserves the equilib- 
rium distribution, the stationary probability distribution 
of the PI-I-GLE scheme requires an accurate integration 
of the stochastic differential equations of motion on the 
timescale of the fastest modes. For this reason, whenever 
a multiple time step integrator is used, the thermostat 
should be applied in the inner loop. As a consequence, 
the integration of the GLE introduces a sizeable compu- 
tational overhead, which can be significant in cases when 
the calculation of physical forces is inexpensive. In the 
case of ah initio molecular dynamics, which is the pri- 
mary target for the present method, the overhead will be 
completely negligible. 



III. THE QUARTIC DOUBLE WELL 



As mentioned in the Introduction, one of the 
fundamental shortcomings of the original quantum 
thermostalP2l is its inability to provide a realistic descrip- 
tion of tunnelling effects. The one-dimensional double 
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FIG. 3. Probability density for a hydrogen atom in a quar- 
tic double well potential with the minima separated by 0.6 A 
and a barrier height of 1000 K. All simulations were performed 
with a target temperature of 300 K. The exact quantum me- 
chanical result (dashed black line) was obtained by numeri- 
cal solution of the Schrodinger equation, the contributions of 
the various eigenstates being averaged with the appropriate 
Boltzmann weight. The four panels compare this exact result 
with the PIMD and Pl-f GLE results with increasing bead 
numbers. 



well potential. 



V{q) = h 



(19) 



in which both zero-point energy and tunnelling are im- 
portant, and can be tuned by adjusting the height li of 
the barrier and the distance d between the minima, there- 
fore provides an ideal first example with which to bench- 
mark the performance of our combined PI+GLE method. 

In Fig. 3 we report the particle density p{q) obtained 
in PIMD and PI+GLE simulations of this potential with 
d — 0.6 A and ti = 1000 K at a temperature of 300 
K, using different numbers of beads. By comparing the 
curves with the exact finite-temperature density 



p(g)=^e-'/'=-^|^,(g)|VE' 



(20) 



where and "0^ are the eigenvalues and eigenfunctions 
of the Schrodinger equation, one sees that the use of a 
GLE dramatically improves the results even when P = 1. 
By increasing P systematic convergence to exact result 
is achieved, and the convergence is greatly accelerated by 
the Langevin dynamics. 

In order to characterise quantitatively the convergence 
of the density, we computed the distance between p{q) 
and p{q), as a function of the parameters of the potential 
and the bead number. For the distance metric we em- 
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ployed the square root of the Jensen-Shannon divergence 



d=0.6 A 



d=2A 



JS 



P) 



p{q) In 



^^'^^^^ p{q) + p{q) ' ''''""Xg) + 

(21) 

which is a metric that is widely used to compare proba- 
bility distributions.^ 

In Fig. 4 we present the convergence of p{q) to p{q) as 
the number of beads is increased. The simulations were 
performed using the mass of a hydrogen atom, different 
values of d and h as indicated in the figure, and a time 
step of 0.1 fs. A very small time step was needed because 
we integrated the PI equations of motion directly using 
the velocity Verlet method, without exploiting the exact 
free ring polymer evolution that becomes possible in the 
normal mode representation.'^l For each set of parameters 
we ran PIMD and PI-GLE trajectories for 20 ns. 

It is clear from Fig. 4 that the PI-GLE simula- 
tions yield a significantly better estimate of the finite- 
temperature density of the quartic double well than the 
PIMD simulations, even in the regime of small d and 
large h where tunnelling plays an important role. The 
addition of the GLE provides a given level of accuracy 
with an effort that is between two and eight times smaller 
than with a standard PI, depending on the parameters 
of the potential and on the accuracy required. 

As the density approaches the exact p{q), all physi- 
cal observables that depend on the position of the par- 
ticle converge to their quantum mechanical expectation 
values. In Fig. 5 we have plotted the average potential 
energy computed from the same trajectories that were 
used to compute the densities in Fig. 4. Again, the use 
of a GLE thermostat consistently improves the quality of 
results. 



IV. LIQUID WATER 

In the previous section we have demonstrated, in 
a simple one-dimensional case, that an appropriately- 
constructed GLE can be used to accelerate the conver- 
gence of path integral molecular dynamics to the exact 
quantum mechanical probability distribution. However, 
as discussed in earlier publications,'^'^ when applying the 
quantum thermostat to a multi-dimensional problem, one 
must be wary of zero point energy leakage. The GLE sets 
the various normal modes to different effective tempera- 
tures, and energy may flow from hot (high frequency) to 
cold (low frequency) modes because of anharmonic cou- 
plings. 

To assess the behavior of the combined PI+GLE strat- 
egy in this respect, and to evaluate the computational 
savings that can be expected from the GLE in a more typ- 
ical application, we have performed some additional sim- 
ulations for liquid water. For these simulations, we used a 
recently-developed flexible water model that was fit to re- 
produce a wide variety of properties of the liquid in path 
integral simulations.'^^ This avoids the double counting of 
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FIG. 4. Distance between the exact density p{q) and the 
densities p{q) obtained from PIMD and PI-I-GLE simulations 
with different bead numbers, for the quartic double well po- 
tential in Eq. (19 1. Six different combinations of parameters 
of the double well potential have been tested. Note that in 
all cases the use of a GLE leads to a significant improvement 
in the estimate of the density. 
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FIG. 5. Average potential energy for a hydrogen atom 
in a quartic double well potential computed by PIMD and 
PI-I-GLE, as a function of the number of beads and for dif- 
ferent parameters of the potential. The exact, quantum me- 
chanical expectation value is marked with a line. 



nuclear quantum effects that would result from the use of 
a force field fit to reproduce experimental data in classical 
dynamics. We performed simulations with both conven- 
tional PIMD and PI+GLE, using different numbers of 
beads. The refined ring polymer contraction scheme of 
Markland et al.-^ was used to reduce the computational 
effort, with the long-range electrostatic interactions be- 
yond 5 A contracted to the ring polymer centroid. The 
use of this scheme has been shown previously to provide 
significant computational savings in empirical force field 
simulations without affecting the accuracy of the results.'^ 
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FIG. 6. The average value of the potential energy and 
the virial kinetic energy for a simulation of a flexible water 
modePD at T = 298 K, plotted as a function of the number 
of beads. The results obtained with conventional PIMD and 
PI-I-GLE are compared, and the value of (V) obtained with 
the original quantum thermostalP^ (QT) is also reported. 



Simulations were performed with a target temperature 
T = 298 K. We used a multiple-time step scheme with an 
outer time step of 0.75 fs and covalently bonded interac- 
tions computed every 0.15 fs. We ran 3.15 ns trajectories 
for each set of parameters, with the first 150 ps used 
for equilibration. Ergodic constant-temperature sam- 
pling was achieved in the PIMD simulations by applying 
a targeted stochastic scheme to the internal modes of the 
necklace^ and a global stochastic velocity rescaling to the 
centroid.!^ 

In Fig. 6 we report the expectation values of the 
potential energy and the centroid virial kinetic en- 
ergy, which were computed using standard path inte- 
gral estimators.^** In this case, the GLE is seen to be 
extremely useful, and leads to much faster convergence 
of averages compared to conventional PIMD. Interest- 
ingly, we found that to converge average energies to an 
error smaller than 1 kJ/mol in the PIMD simulations, 
P should be increased well beyond the 32 beads that 
are generally adopted for room-temperature water. Con- 
versely, PI-I-GLE yields results in perfect agreement with 
a 128 bead PIMD simulation using fewer than 16 beads. 

As a more sensitive benchmark of the convergence 
of the properties of quantum water we have also com- 
puted the constant- volume specific heat cy, which is 
known to require very large values of P for an accurate 



determination.'22l The value of and statistical uncertainty 
in cy were obtained from a quadratic fit to the values 
of {V) and (T) computed at 293, 298 and 303 K. The 
PI+GLE method is again seen to perform exceedingly 
well in this test (see Fig. 6). The convergence is some- 
what accelerated by a fortuitous cancellation between the 
errors in the Pl-f GLE estimates of the potential and ki- 
netic energy contributions to cy , which are however both 
very well converged by the time P ~ 8. 

Comprehensive results for the convergence of struc- 
tural properties of water are reported in Fig. 7, in which 
we systematically examine the radial distribution func- 
tions that are obtained with different bead numbers. 
Here again the GLE helps tremendously in accelerating 
the convergence of the PI calculation, in particular in the 
region of the intramolecular peaks, which are strongly af- 
fected by nuclear quantum effects. The only case in which 
the Pl-f GLE simulation is in worse agreement agreement 
with the exact result than the PIMD simulation is when 
the oxygen-oxygen g(r) is computed with too few beads 
{P — 1 and P = 2). This is a consequence of the zero- 
point energy leakage in the PI-I-GLE simulation, which 
heats up low frequency degrees of freedom and washes out 
the long-range features from 5oo('')- By increasing the 
strength of the thermostat in the low-frequency region 
(as is done for instance in the case of the original quan- 
tum thermostat - see Sec. II. D) it is possible to mitigate 
this effect, at the expense however of a greater distur- 
bance on diffusive modes, and hence on the efficiency of 
sampling. The PI+GLE strategy allows one to counter 
the effect of zero point energy leakage with a moderate 
increase in the number of beads, and indeed the PI+GLE 
radial distribution functions for liquid water are seen to 
be significantly more accurate than the PIMD radial dis- 
tribution functions in Fig. 7 for all P > 4. 

Combining the results in Figs. 6 and 7, one sees that 
PI+GLE provides the same accuracy in the thermody- 
namic and structural properties of liquid water with 8 
beads as PIMD provides with 32 beads. This is the 
level of accuracy that is commonly accepted for room- 
temperature water simulations, and adding an appropri- 
ate GLE to the PI simulation allows it to be achieved 
with a four-fold reduction in the number of beads. The 
comparison becomes even more favourable if higher ac- 
curacy is required, with a 16-bead PI+GLE simulation 
being as accurate as a 128-bead PIMD simulation for all 
of the observables we have considered. 



V. CONCLUSIONS 

In the present paper we have discussed and thor- 
oughly demonstrated how a properly-designed general- 
ized Langevin equation can be used to accelerate the 
convergence of path integral molecular dynamics to the 
exact quantum mechanical thermal expectation values. 
This leads to substantial savings in computational effort, 
the number of beads required to obtain a given level of 
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FIG. 7. Errors in the radial distribution functions of liquid water at T = 298 K obtained in low-P simulations using PIMD and 
PI+GLE, relative to a fully-converged PIMD reference simulation with 128 beads. Note that the scale of the y-axis is greatly 
magnified as P is increased and the Ajr's become smaller. The results obtained using the original quantum thermostalP^ (QT) 
are also reported in the panels with P = 1. The statistical error in the distribution functions is of the order of 10"'^. 



accuracy being reduced by a factor of 4 or more. The 
original (one-bead) quantum thermostalP^ provides an 
even cheaper way to include nuclear quantum effects in 
mildly anharmonic systems, but it is an inherently ap- 
proximate technique. The present PI-I-GLE combination 
captures tunnelling effects as well as zero point energy 
effects, and it can be systematically improved to give the 
exact quantum mechanical result simply by increasing 
the number of path integral beads. We expect that this 
combination will be particularly valuable when used in 
conjunction with ab initio molecular dynamics, in which 
the forces acting on the nuclei are so expensive to evalu- 
ate that nuclear quantum effects have only seldom been 



considered in the past. 

One final observation is that we have concentrated ex- 
clusively in this paper on the standard second order Trot- 
ter product path integral in Eqs. Q and (4]). Another 
way to reduce the number of path integral beads is to 
use a higher-order imaginary time propagator, and some 
interesting progress has been made in this direction over 
the years. ■^^ '^^ This approach is fundamentally different 
from the approach we have investigated here, and poten- 
tially complementary to it. One could in principle de- 
velop a GLE scheme to accelerate the convergence of any 
imaginary time propagator, including the more promising 
of the higher-order propagators that have been suggested 
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in the literature. In this way, one might hope to be able 
to reduce the number of beads even further, and make 
the inclusion of nuclear quantum effects in ab initio simu- 
lations almost routine. In any event, the results we have 
presented in this paper clearly demonstrate the potential 
of the generalized Langevin equation as a computational 
tool for accelerating the convergence of PIMD. 



Eq. ( 13 ) and the non-uniqueness of its solutions, and Gio- 
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